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=2 ! Abstract 

B 

pH ■ We have employed the method of spectral moments to study the density of 

o ■ 

O ■ vibrational states and the Raman coupling coefficient of large 2- and 3- dimen- 



sional percolators at threshold and at higher concentration. We first discuss 



■ the over-and under-flow problems of the procedure which arise when -like in 

the present case- it is necessary to calculate a few thousand moments. Then 
we report on the numerical results; these show that different scattering mech- 
anisms, all a priori equally probable in real systems, produce largely different 
coupling coefficients with different frequency dependence. Our results are 
compared with existing scaling theories of Raman scattering. The situation 
that emerges is complex; on the one hand, there is indication that the existing 
theory is not satisfactory; on the other hand, the simulations above threshold 
show that in this case the coupling coefficients have very little resemblance. 
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if any, with the same quantities at threshold. 
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I. INTRODUCTION 



Crystalline solids are characterised by translational invariance and this symmetry prop- 
erty makes it possible to derive explicit analytic expressions for the vibrational eigenvectors 
in the harmonic approximation or, at worst, entails the diagonalization of small matrices. 
As a consequence, a lot is known on many physical properties that are determined by the 
vibrational dynamics of crystals (vibrational density of states, thermal properties, inelastic 
light and neutron scattering, and so on). The situation is completely different in disordered 
solids: in this case the lack of translational invariance often requires that the dynamical 
properties be studied by numerical calculations on more-or-less realistic models. 

Even when the model well represents the real system, one of the factors that make a 
numerical calculation useful is the dimension of the sample; in many cases, the bigger the 
dimension, the lower the probabihty that finite size effects dominate the results and make 
any comparison with experimental data problematic. 

As regards the study of the vibrational dynamics of disordered systems, and assuming 
that the harmonic approximation holds, one possible numerical approach is to build up the 
dynamical matrix and diagonalize it; this provides vibrational eigenvalues and eigenvectors, 
i.e. all the information that may be required. The problem is that the linear size of the 
matrix to be diagonalized is as large as the number of vibrational degrees of freedom, and, 
for example, a sample containing 10x10x10 atoms would require the diagonalization of a 
3000 X 3000 matrix: a very large matrix for not too large a system. Even though in some 
cases sufficient insight may be gained by considering only one degree of freedom per atom, 
the calculations are in general lengthy and expensive. In any case, 3-dimensional models 
containing, say, 20x20x20 atoms are at present practically intractable with this technique. 

In many cases all the detailed information contained in the eigenvalues and eigenvectors 
is not necessary. For example, the density of vibrational states is a more significant quantity 
to compare with experiment than the sequence of eigenvalues. The important thing is that 
some of these relevant quantities can be calculated without diagonalizing the dynamical 
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matrix, by exploiting in a different and more efficient way the information it contains. 

In the present paper we will use the method of moments to compute the density of 
vibrational states and the Raman coupling coefficient of model disordered systems consisting 
of two- and three-dimensional (hereafter 2D and 3D respectively) site and bond percolators; 
these were studied both at percolation threshold concentration and at higher concentration. 
Percolators at threshold are fractals and using arguments based on scale invariance it has 
been shown that the density of vibrational states follows a power law of the type 
p{u) oc u!'^~^, where d is a parameter known as the spectral dimension. 

In the recent literature several attempts were made to find a similar power law for the 
Raman coupling coefficient C{uj), which determines the (Stokes) scattered intensity in the 
following way: 

/(cu) = (n(cu) + l)C{u;)p{u;)/iJ (1) 

where n{uj) is the Bose-Einstein population factor. As we will discuss in the following 
sections, the proposed scaling laws C{uj) oc are based on assumptions that are not 
universally accepted so that comparison with numerical simulation is of great importance. 
Numerical calculations of C (uj) were earlier produced for a number of percolators at threshold 
10-0] using dynamical matrix diagonalization. As mentioned, this allows for only limited 
dimensions of the models, and this restriction is particularly bad in this case, since there are 
indications |0,H that power laws for C {uj) might apply only at low frequency. In the present 
paper we present an exhaustive analysis of C{uj) in much larger systems, in order to avoid 
as far as possible finite size effects and reach the lowest possible frequency. 

In section II we summarize the method of calculation; section III is devoted to the 
stability and precision problems of the procedure and to the way we were able to overcome 
over- and under-flow difficulties; in section IV we present the numerical data on various 
types of percolators; the data are discussed in connection with the phenomenology of Raman 
scattering from real disordered solids in section V. 
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II. THE METHOD OF MOMENTS 



In the following we shall outline the method of spectral moments, which is reported in 
detail in the papers by Benoit and coworkers |]T[]. This method applies to the case where one 
is interested in calculating a spectral function /(cu^) of the form: 

3N 

f{u^') = Y.dl6{uj'-ujl) (2) 

A=l 

where 

3N 

d\ = ^Piei{X) (3) 

i=l 

The system consists of N identical masses which interact via identical harmonic potentials 
and is characterized by the 3A^ x 3A^ dynamical matrix D with eigenvalues uj\ and eigen- 
vectors ej(A); A = 1...3A^ labels the normal modes and z is a collective index that labels the 
masses (/ = l...A^) and the cartesian components (a = x, y, z) of the eigenvectors, i = {I, a). 
The coefficients Pi depend on the spectral function to be computed; the explicit expressions 
of Pi for the cases of interest in this paper are discussed in the Appendix. 

It has been shown by Benoit and coworkers that f{oj'^) can be put in the form 

fi^") = -l limIm{R(z)} (4) 

where z = oj"^ + is, and 

J-OD Z — OJ^ 

The latter can be developed as a continuous fraction 

R{z) = ^—^^ (5) 

z-ai r 



z - a2 

z-as- 

where the real coefficients a„ and 6„ depend on the generalized moments Unm and Vnm of 
f{uj'^) (whence the name of the method): 



where 



nn 



and 



nn 



C P^{oo')Pr,{uj^)f{oo')du;' (7) 

JO 



I' P^{uj')Pr,{uj^)f{uj')uj^ du;\ (8) 

JO 



{Pn{oj'^)} is the succession of polynomials orthogonal with respect to the spectral function 
f{uj'^) (as in the rest of this paper, the frequency scale in equations 7 and 8 has been 
normalized so that the maximum frequency is 1). The moments in turn can be computed 
by a recursive procedure that makes use of the dynamical matrix. Let t° be the normalized 
vector with components proportional to Pi, it turns out that 

z/„„ = (t(-),t(")), I7„„ = (t(-),Dt(")) (9) 

where t^"-* is a vector which obeys the recursive relation 

= (D - a„+i)t(") - 6„t("-i) (10) 



Equations || to |10], with t*^"^) = 0, determine a„ and 6„ recursively, and these give /(cu^) 
through R{z). 

In principle, if one computed a number of a„ and 6„ coefficients equal to the number of 
distinct eigenvalues |^ , in the limit e ^ one would reproduce exactly the response function 
in Eq. (§. 

In computations of practical interest (i.e. on systems with a large number of masses) this 
is neither possible nor necessary. The number of moments to be computed and the value 
of e depend on the frequency range where the spectral function is to be reproduced most 
carefully. 

The use a finite value for e, which is equivalent to considering the convolution of /(cu^) 
with a lorentzian 



makes it impossible to resolve spectral features with width of the order of (or smaller than) 
e. 

The fact that only a limited number of moments is computed introduces a truncation 
error that, as noted by Turchi et al. |0, can be minimized when the spectral function /(cu^) 
has no gaps and no divergences. In this case the coefficients a„ and 6„ tend to make small 
fluctuations about constant asymptotic values as n oo, and this suggests to substitute 
the asymptotic values in Eq. (|^) for the coefficients that are not computed. If one wants 
only M moments one re-writes Eq. as: 
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R{z) 



hi 

z — ai 



z - a2 - ■ ■ 



z-ttM- Tm{z) 

and under the assumption of constant an=a and for n > M, Tm{z) is easily found to 
obey 



z — a — hTM{z) 

Since in our case a„ and 6„ continue to make small oscillations about a and h respectively for 
large n, this procedure provides only an approximation to the true R{z)] in any case Tm{z) 
becomes less and less important as the number of moments that are actually computed 
increases. 



III. STABILITY AND PRECISION OF THE PROCEDURE 

Since we will be mainly interested in studying the behavior of spectral functions like 
in Eq. (|^) at frequencies of the order of 10^'^ ^ 10~^ times the Debye frequency, the use 
of many (some thousand) moments is necessary, as can be seen empirically and deduced 



theoretically. This requires many iterations of Eq. ([ToD that can introduce numerical errors 
in the computed /(cu^). The errors can arise both from the precision of the calculation 
(accumulation of roundoff approximations) and/or from the instability of the algorithm, i.e. 
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from the sensitivity of the algorithm to small random departures of and 6„ from the true 
values. 

As a first example in Fig. 1 we show the effect of the number of computed moments on 
the density of states of 2D, 650x650 site percolators at percolation threshold, consisting of 
identical masses connected by identical springs; each mass has only one degree of freedom, 
cyclical boundary conditions are imposed [[l^. It is clearly seen that in passing from 50 to 
5000 iterations the minimum frequency at which the expected density of states {p{uj) oc u^'^^) 



Tl| is reproduced, decreases by more than one order of magnitude. Note that 1000 iterations 
are not sufficient. The origin of the dip at uj ^ 0.7 will be discussed later. 

The necessity of so many iterations may be justified theoretically as follows. Using the 



orthogonality relation 







with k = 0, ...,n — 1, equation (0) becomes: 



iynn= / Pn{uj')iuj'-l/2rf{uj')du;' 







from which it is easily seen that, for large n, Vnn is mostly determined by the extremes of the 
uj range, so that the new information that is added at each step is more and more concerned 
with the extremes; the new information is necessary to approach the minimum attainable 



frequency, as seen in Fig. 1 |T^. Moreover, it is also clear that in order to reproduce the 
behaviour of fiuj"^) in the central region of frequencies it is neither necessary nor (as we will 
see) advisable to use many iterations. 

The first numerical difficulty that is met in the calculation of many moments is that z/„„ 
(and Vnn) tend to approach or oo, causing overfiow or underfiow of the computed moments. 
This difficulty can be overcome by noting [||] that if one scales the physical characteristics 
of the system (masses and/or elastic constants) by a factor s, the generalized moments scale 
as well: 
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a„ ^ s^a„ (11) 

It is therefore possible to choose a value of s such that the generalized moments are finite. 
Since 6„ = z/„„/ z/„_i one possibility is to start with an arbitrary elastic constant /Ci, then 
compute a sequence of M\ coefficients 6„ until vm^Mx reaches the over /under- flow condition, 
then compute (6) = Mf ^ Z^n^n and scale /Ci in agreement with the first of equations |ll| so 
that (6) is equal to a desired value. 

At first sight, it would seem that (6) = 1 is a good choice to avoid divergences, so that 
for the scaling factor s we have 

s^ = ^ 
ib) 

This procedure may be repeated several times, until the desired number of moments can be 
calculated. 

However (b) = 1 is not the most convenient choice. In fact, since the 6„'s fiuctuate 
around (6), we can put 6„ = + with (e^) = 0; in this way we have 

or 



n 
i=l 

= (6)Voo(l + E^^ + E^^^. + -) (12) 

i i^j 

Since Z]j = 0, we have 

ijtj ij i i 

SO that 

i^nn ~ ^^oo(fe)"(l - n{e^)) ^ z/oo(fe)"e-"<^'> (13) 



Thus, in order to control the divergence of z/„„ one should not impose (6) = 1, but 

(b) = e^^'l (14) 

We have verified that with this choice the number of corrections of the elastic constant 
necessary for all the 5001 iterations of a medium-sized system (square lattice of linear size 
100), or for 5000 iterations of a larger system (2D percolator of linear size 650), is smaller 
than with the first choice. 

However, neither correction procedure eliminates a further problem concerning the or- 
thogonality of the series {t^**)}. In fact, from Eq. (p!0|) it is easy to see that (t^'^\t^"^^) oc 5nm, 
but in practical calculations the orthogonality is lost if \n — m\ is large enough. 

In Fig. 2 we report (t^^^t^"^) as a function of n for four different systems: full (ordered) 
square lattices of linear sizes 50 and 100, and site percolators at threshold concentration of 
linear sizes 50 and 650. The calculations were carried out in double precision, but using 
single precision only increases the value of the initial plateau. It is interesting to note that 
the value of n at which orthogonality begins to fade does not depend in a simple way on the 
number of masses (and therefore on the number of required computer operations): in fact, 
for the larger ordered lattice orthogonality is maintained for a greater number of iterations, 
while for percolators this effect is much less evident, if any. 

This behaviour can be understood by noting that the recursive relation ^ is very similar 
to the one used in the Lanczos procedure of tridiagonalization, and in fact the two become 
identical if in the Lanczos procedure [|1^] one makes the substitutions: 



a — >■ a 
(3^^b 



As discussed at length in reference |]T3[ , this loss of orthogonality is inherent and not caused 
by accumulation of numerical roundoff errors. Moreover, it is known that the Lanczos pro- 
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cedure without any re-orthogonalization works well for the extreme eigenvalues and eigen- 
vectors, while there are problems for the inner range. 

From another point of view, the iterative equation (|1^) may be considered as a non linear 
logistic application (the non linearity coming from the dependence of a„ and 6„ on t*^")) in 
an iV-dimensional space. It is known that these equations admit chaotic solutions for the 
t(") variables, i.e. they produce sequences {tW} that, as n grows, become very unstable 
with respect to little variations of t'-''-*. Therefore, for example, from t'-'^-* numerically we 
obtain t''-^-* = t'-^-* + ^t*^^-*, ^t'-^-* being the machine precision error. Consequently, because of 
the intrinsic chaoticity of equation (p!0|) , the computed {t'^"^} sequence will diverge from the 



"true" one and (t(°),t'(")) ^ 0. 

Our numerical data are in complete agreement with the above results of matrix calculus. 
First of all, when we compute the density of states (or the Raman coupling coefficient) of 
a linear chain, orthogonality is preserved no matter how long the chain; this is expected 
because in this case the dynamical matrix is tridiagonal. 

Moreover, when comparing the outcome of the method of moments with the eigenvalues 
of an exactly solvable model (square harmonic lattice with identical masses and springs with 
fixed boundary conditions |T^) we find that as orthogonality begins to be lost the peaks of 
the density of states no longer fit the exact eigenvalues in the central part of the (linear) 
spectrum, while the extreme eigenvalues are still well fitted, provided we use a sufficient 
number of moments. In any case, the central part of the spectrum is perfectly fitted only for 
systems small enough that computation of all moments does not result in non-orthogonal 
vectors. For larger systems, increasing the number of moments worsens the look of the 
central part of the spectrum, be it the density of states or the Raman coupling coefficient, 
and this is the origin of the dip in the density of states of Fig. 1. In our case, this is not a 
serious problem. In fact, on the one hand we are mostly interested in the low frequency part 
of the coupling coefficients, this being the very reason why we compute so many moments; 
on the other hand, if we should be interested in the central range we would compute only a 
few moments. This is possible because the coupling coefficients are smooth functions due to 
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disorder, we don't need to resolve one peak from the next, all we need is a frequency- averaged 
result. 



IV. NUMERICAL RESULTS ON RAMAN COUPLING COEFFICIENTS 

We proceed now to presenting the numerical data. As mentioned in the introduction, the 
problem of finding a scaling law for the Raman coupling coefficient in fractals, if it exists, 
has been the object of several experimental and theoretical papers ||7|,p!5|- [l9|] . Besides being 



interesting on their own, fractals are thought to be reasonably representative and relatively 
simple models of important classes of disordered systems, especially as regards vibrational 
characteristics. Raman scattering is one of the most convenient experimental techniques to 
study vibrations in disordered solids: it does not require big facilities and, contrary to the 
case of crystals, all vibrational modes contribute to the scattering in a disordered system. 
The price to be paid is that C{uj) in equation (|^) is not known a priori, so that, as mentioned, 
an important question regards whether or not it can be cast in the form C{uj) oc uj^, and 
what expression x has in terms of the spectral dimension d, of the fractal dimension D, and 
possibly of other parameters which characterize the fractal. 

All the expressions proposed for x so far were derived under the assumption that C (u) can 
be determined by considering the scaling properties of the strain induced by the vibrations 
of the fractal (the so-called fractons) [|^,|T5|,|16|,[18|-^ . In two recent papers we argued 



against such possibility on the basis of numerical calculations of C{u!) in 2D and 3D site 
percolators, which apparently showed that (z) none of the proposed x could fit the numerical 
data in 2D, and (ii) a single x seems not to be sufficient in 3D for the Dipole- Induced-Dipole 
(DID) scattering mechanism (see Appendix). 

The conclusion we arrived at on the basis of these numerical results, i.e. that the existence 
of scaling for C{uj) is not always evident in fractals, has been recently challenged on the basis 
of the following arguments: (z) the size of the percolation clusters employed (60x60 in 2D, 
29x29x29 in 3D) was too small to allow positive conclusions to be drawn on scaling, which 
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is expected to hold at low frequency and (ii) even more so for s2te-percolators: bond- 
percolators should be used instead, because for the latter the scaling regime starts at higher 
frequency 0. 

These issues will be discussed in the following, after we have reported numerical data 
obtained with the method of moments for a variety of much larger systems including site- and 
bond-percolators at and above percolation threshold, whose scattering units may or may not 
have random electrical polarizabilities (electrical disorder), and assuming different scattering 
mechanisms: DID, DID truncated to nearest neighbors (NNDID), and Bond-Polarizability 
(BPOL). 

In Fig. 3 we report C{uj) for the DID and BPOL scattering mechanisms for 650x650 
2D site percolators at threshold. The interesting features in this figure are: (i) The BPOL 
spectrum shows a single slope, m = 1.24, down to the minimum frequency, (ii) the DID 
spectrum clearly shows a crossover at ~ 8 x 10~^ from m ^ 0.94 (high frequency) to 
m ^ 0.8, which was uncovered in due to limited size and statistics. 

In Fig. 4 we show the same coupling coefficients for 3D site percolators of linear size 
80. In agreement with the results of diagonalization, P,|25| there is a very evident crossover 



in the DID case, while the BPOL spectrum again is a straight line in almost the whole 
frequency range. The open circles represent the DID spectrum of clusters of linear size 
40. Apart from the obvious fact that the minimum frequency is greater in this case, we do 
not observe any significant difference with respect to the larger clusters, indicating that, at 
least for C{uj), finite size effects play a minor role, if any. It is interesting to mention that 



removal of dangling bonds ||2^ in the site percolators does not change the look and the slope 
of the BPOL coupling coefficients both in 2D and 3D, but it does change the density of 
states which does not appear to follow an uj^ law in the whole frequency range; this effect 
is particularly evident in 3D. 

The C {uj) 's relative to 2D (linear dimension 500) and 3D (linear dimension 70) |2^ bond 
percolators are reported in Figs. 5 and 6 respectively, for both scattering mechanisms DID 
and BPOL. The BPOL spectra are well fitted by straight lines with almost the same slopes 
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(m2D = 1-26, m^D = 1-58) as for site percolators. The DID spectra have a roundish look, 
but in the low frequency part they are reasonably fitted by straight lines. 

For bond percolators, The nearest-neighbor DID coupling coefficients in 2D and 3D are 
practically identical to DID, in agreement with what found in ref. and we will not report 
them. This fact can be easily understood considering that in bond percolators (contrary to 
site percolators) it is possible to find nearest neighbor sites which are not connected by any 
bond. This results in highly uncorrelated motions of nearby masses that produce most of 
the DID scattering. 

In Fig. 7 we show the effect of the so called electrical disorder on the DID scattering 
of 2D site percolators. Each mass is randomly assigned one of two values ai or a2 of the 
bare polarizability; the curve of Fig. 3 (ai = ^2 = 1) is also reported for comparison. Note 
that spectra with different electrical disorders coalesce at low frequency {u < 2 x 10^^). 
Qualitatively the same result is obtained with all other types of clusters, with and without 
dangling bonds. In the case of BPOL for 2D site percolators the C(a;)'s with electrical 
disorder show a change of slope at u ^ 10~^. 

In Fig. 8 are reported the DID and BPOL spectra for 3D site percolators having a 
concentration (c = 0.5) higher than the percolation threshold concentration, together with 
the relative density of states. The latter exhibits the well known |^ crossover at Up ^ 
7 X 10~^ with change of slope. The fact that neither slope is equal to what expected for 
phonons (m = 2) or fractons (m = 0.31) is also known [^. The C{uj)'s also change their 



frequency dependence at a frequency, uc ~ 6 x 10^^, which is apparently slightly lower than 
Up. The spectra are rather noisy and it is impossible to judge whether C{uj) is a straight 
line at low frequency, but in any case it is clear that below uc the "slopes" are very different 
from what was found for the same percolators at threshold, see Fig. 4. The same qualitative 
behavior is observed in 2D, and for bond percolators, though the 3D spectra in this case are 
even more noisy at low frequency. 

The slopes observed are summarized in Table I for the different kinds of percolators and 
scattering mechanisms considered; in the case of DID, where different slopes ore found in 

14 



the same spectrum, we report the low frequency ones for the reason discussed in the next 
section. 



V. DISCUSSION 

From the results presented in the previous section, it appears that C{uj) behaves quite 
differently depending on the scattering mechanism. In fact, in the case of BPOL the nu- 
merically computed coefficients follow very closely the law C{uj) oc uj^ almost in the whole 
investigated frequency range; the values of x do not depend appreciably on site- or bond- 
percolation, nor on the presence or absence of dangling bonds. The situation is rather dif- 
ferent for DID scattering; with this mechanism, in none of the case studied does C{uj) oc uj^ 
hold in the whole range of frequencies. Following the arguments of refs. we shall as- 

sume that it is the low frequency part of the DID spectra that is to be compared with the 
proposed scaling laws (but it is not clear to us why the arguments of those references do 
not apply to BPOL). This is supported by the result of Fig. 7, that shows that the effect of 
electrical disorder, which is a local random disturbance and as such is not expected to scale, 
disappears at more-or-less the same frequency {u ~ 10^^) where C{uj) changes slope. 

Since the DID and BPOL coupling coefficients are so different and produce different 



slopes, we can compare our results only with the theoretical model of Alexander et al ||T9 
because to our knowledge this is the only paper where different scattering mechanisms are 
considered. Alexander et al consider DID and NNDID: the latter coincides with BPOL for 
site percolators and, as mentioned, it is practically identical to DID for bond percolators, 
so that it is possible to compare our numerical data with the formulas of ref. [|19i that are: 



Ciuj)DID OC ,;(2'i/^)(-+'^)-3'^ (15) 

and 

C{u)mndid oc u;^^'^/^ (16) 
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where a is a scaling index introduced by Alexander et al whose value is not known 



priori: for a homogeneous medium cr = 1, while a < \ means a violation of scaling [jl9 
Values of 0" ~ 1 were deduced in ref. on the basis of simulations 0] performed on smaller 
clusters than done here. Equation (^6|) was first derived by Boukenter et al |15[ . 

Fitting the slopes with equations (|T^) and (|TB|) we obtain the values of a reported in 
Table I. The picture that emerges is complex. 

• {i) In two cases, 3D site percolators (DID) and 2D bond percolators (DID) we find 
cr ~ 1, which is the desired value; 

• {ii) in two cases, 2D site percolators (NNDID) and 3D bond percolators (DID), the 
values found are smaller than 1, but not too much (a 0.9); 

• {Hi) in the remaining four cases, a is either too large to be plausible (a ~ 1.4 for 2D 
site percolators (DID) and 3D site percolators (NNDID)) or exceedingly small (cr 0.3 
for 2D and 3D bond percolators (NNDID)). 

The first thing that we note is that so far a has always been thought of as a parameter 
that characterizes the scattering system, not the scattering mechanism, so that the wide 
variations observed in Table I between DID and NNDID apparently imply that equations 



(15) and/or ([Tq) are not right. 

On the other hand, the fact that the BPOL C(co')'s are straight lines over a frequency 
range of more than two orders of magnitude, and are so insensitive to the nature of the 
system, might reasonably be taken as indication that for this scattering mechanism scaling 
holds, but following a law other than equation ([T6|) for site percolators. In our opinion it 
would not be too much of a surprise if BPOL should scale and DID should not because 
the modulation of polarization in DID proceeds through electromagnetism, which definitely 
does not propagate along the fractal paths. In any case, we think that nothing should be 
taken for granted when treating this subject. For example, as mentioned in the previous 
section, removal of dangling bonds from 2D and (especially) 3D site percolators produces 
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a situation where the mass distribution follows the same power law as in the case with 
dangling bonds, Cbpol{^) does the same, but the density of states is no longer a power 
law for uj >^ 4 x 10^^. This situation might even be interpreted as indication that BPOL's 
being a straight line has actually nothing to do with dynamic scaling: in fact, in these 
systems Cbpol{^) is a straight line in a frequency region where the density of states (i.e. 
the most believed scaling dynamical quantity) is not. At present we are not able to resolve 
this question. 

Another result in Table I that is worth a comment is the value of a for 2D and 3D bond 
percolators under the NNDID scattering mechanism. The fact that for these percolators 
DID and NNDID C{uj)'s are almost identical implies that at all frequencies available to the 
present simulation the scattering coefficient is mostly determined by pairs of nearest neighbor 
masses that are not connected by a bond. This situation produces cr ^ 0.3 for NNDID. In 
principle, however, it cannot be excluded that on much larger clusters the contribution of 



these pairs may become less important. In order to find a ^ 1 the simulation should 
yield NNDID slopes of ~ 1.1 (3D) and 1.4 (2D), i.e. very different from the values obtained 
here. In any case, these new slopes would be found in a frequency range (cu ~< 10""^) so 
low to be of little physical significance. 

This remark introduces another important issue, i.e. the possibility of extracting from 
experimental spectra information on the static and dynamic fractal parameters. This is in 
general done by fitting the low-frequency parts of the spectra (where scaling is expected to 
hold) with expressions like Eqs. (plSj) or (p^G]). Assuming that one has the right equations, 
one obvious shortcoming of this procedure is that one should know the scattering mecha- 
nism (DID, BPOL, NNDID, a mixture, ...), which has been shown to be so important in 
determining the low frequency slope. But even if this was known, there is another problem 
that is evident from fig. 8. In fact, real disordered materials cannot reasonably be thought 
of as percolators at threshold: if the percolation model is to have a sense, it must be a 
percolator above threshold. From fig. 8 we see that in this case the values of the "slopes" 
bear no resemblance with those found at threshold, and especially so at low frequency. This 
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fact has long been known as regards the density of states of percolators above threshold p6 



Of course, the particular concentration value used here (c=0.5) has no special meaning, nor 
is it thought to represent any particular physical system. 

In conclusion, we think that the search for the scaling laws which might possibly gov- 
ern Raman scattering from fractal objects is interesting in itself and as such worth being 
pursued; on the basis of the present numerical simulations, we think that these laws have 
not been found yet. Even more disputable is the extraction of fractal parameters from the 
experimental Raman spectra. In this regard, we frankly hope that the present paper may 
encourage a realistic reconsideration of previous work, even if this might lower the role of 
the fractal model of disordered solids. 
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APPENDIX: DEFINITION OF THE INITIAL VECTORS 

The general expression for the Raman intensity in a, (3 polarization at energy hu and 
exchanged momentum is given by: 

I^,{co) =Af dt e^-'Y.{<P^)<M e^''-(^'«-^"(°))) (Al) 

IV 

where ^ is a constant and vr^^(t) is the a, j3 component of the effective polarizability tensor 
of the Z-th atom placed at site R/(t). The time dependence of this tensor is caused by the 
relative displacement of atom / with respect to all other atoms. 
The instantaneous position Ri(t) may be written as 

Ri(t) = + u,(t) 
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where u/(t) is the displacement from the equihbrium position x;. For small displacements, 
both the effective polarizability and the exponential function in equation (|A1| ) can be ex- 
panded in power series: 

gik.R,(i) ^ g^k-x, (1 + ■ U;(t)) (A2) 

vri,(t)-^i, + EE ^ W^)- 
The displacement, in turn, may be expanded in normal modes: 



um^it) = ^^^T. ^e„^(A)A(t) (A3) 



where ^a(^) are the normal coordinates that (for the Stokes part of the spectrum) obey the 
relation: 

/ dt e'-\A,{t)Ay{0)) = + _ (A4) 

J UJ 

From the previous equation we have 

/(^) = ^[!!M±ilp(^)c(cu) (A5) 



having defined 



C{uj) = J2 \Cx\'S{uj - u^)/ J2 ^(^ - ^a) (A6) 



Y,\CxMu-ujx)/p{^ 

A 



and 



Ca = ^EE|^KW-en..(A)]. (A7) 

In writing the above equations we have assumed that we are interested in the depolarized 
component of the spectrum, (a,/3) = {x,y). Considering that 

E^ = o 

equation (|A7|) may be cast in the form 
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or 



= ^EE^e.7(A). (A8) 



m 7 



where 



C\ in equation (^) has the same structure as dx in equation (H), so that the components 
of the initial vector t^^^ are Pi = Pm-y- We note that from equation (|A6|) we obtain 

Pj^^ = E \c.?'-^^ = E ic.r^*(-^ - -I) (All) 

A A 

to be compared with equation (0) in the text. 

In this paper we have considered different mechanisms of polarizabihty modulation, i.e. 
different expressions for vr^^, that result in different t^°^'s. In particular, the explicit expres- 
sions for the studied scattering mechanisms are: 

DID: nlf, = ^Tg(/m)«,«™ 
NNDID : Trlf, = Tg(/m)aza™ 

m,{l} 

BPOL : vri^ = E TS(/m)V^„ 

m,{l} 

where ai is the bare polarizabihty of atom I, T^^(r) is the dipole propagator 

T£U..«„(r) = ViV2...V.^ 

the index m, {1} indicates that the summation concerns all the m nearest neighbors of atom 
/, and Vim is 1 or according to whether atoms / e m are connected by a bond. 
Therefore, from equation ( [A1U| ) we find that: 



DID: = ETS,(/m) 



aiar 
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NNDID : p, = Tg^(/m)«z«^ 

l,{m} 

BPOL : p, = Y: T(5^(/^)^n. 

l,{m} 

In the case of DID the effective polarizability depends on the positions of all atoms, each 
weighted by the factor l/R^; in the NNDID case this effect is limited to nearest neighbors, 
to simulate induction mechanisms that decay faster than 1/ E? . In the case of BPOL mod- 
ulation of polarizability occurs only if the nearest neighbors are actually connected by a 
bond. Therefore, in site percolators NNDID and BPOL coincide, while they do not in bond 
percolators. 

We have also been interested in the calculation of the density of states; in this case the 
values of Pi are uniformly and randomly distributed between -0.5 and 0.5, as shown in ref. 

i- 
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FIGURES 

FIG. 1. Log-log density of states of 2D site percolators (average of 10 realizations) having linear 
dimension L=650, with e = 2 • 10~^ and: (a) 5000 moments; (b) 1000 moments; (c) 50 moments. 
Lo = Lo/umax', iTi is the fitted slope; the curves are shifted vertically for graphical convenience. 

FIG. 2. Scalar product (t(°),t(")) (see text) as a function of the number of iterations n for 
ordered square lattices of linear sizes 50 (a) and 100 (b), and 2D site percolators at threshold of 
linear sizes 650 (c) and 50 (d). 

FIG. 3. C{uj) for 2D site percolators, L=650. (a) DID; (b) BPOL. Average of 10 realizations. 

FIG. 4. C(a;) for 3D site percolators, L=80. (a) DID; (b) BPOL. Average of 20 realizations. 
Open circles: DID for 3D site percolators, L=40, average of 100 realizations. 

FIG. 5. C{uj) for 2D bond percolators, L=500. (a) DID; (b) BPOL. Average of 10 realizations. 

FIG. 6. C{uj) for 3D bond percolators, L=70. (a) DID; (b) BPOL. Average of 10 realizations. 

FIG. 7. Effect of electrical disorder on C{uj) for 2D site percolators, L=650. (a) ai=0, a2=2; 
(b) ai=0.5, a2=1.5; (c) ai=a2=l. Average of 10 realizations. 

FIG. 8. C{uj) for 3D site percolators, L=70, above percolation threshold (c=0.5). (a) Density 

of states; (b) DID; (c) BPOL; the curves are vertically shifted for graphical convenience. Average 
of 10 realizations 
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TABLES 

TABLE I. Fitted slopes, m, for C{uj), see Figures (3) to (8), and values of a obtained by fitting 
the numerical data with equations (15) and (16) for DID and NNDID respectively. For BPOL no 
such theoretical expression is available, but in the case of site percolators BPOL coincides with 



NNDID. 







DID 




NNDID 


BPOL 






m 


a 


m 


a 


m 


site 


2D 


0.80 ±0.02 


1.40 ± 0.01 


1.24 ±0.01 


0.89 ± 0.01 


1.24 ±0.01 


percolator 


3D 


0.41 ± 0.03 


0.97 ±0.03 


1.56 ±0.01 


1.40 ± 0.01 


1.56 ±0.01 


bond 


2D 


0.40 ± 0.05 


1.12 ±0.04 


0.40 ± 0.05 


0.30 ±0.05 


1.26 ±0.01 


percolator 


3D 


0.30 ± 0.08 


0.86 ± 0.07 


0.30 ± 0.08 


0.30 ± 0.05 


1.58 ±0.01 
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